DISTRIBUTION  STATEMENT  A.  Approved  for  public  release;  distribution  is  unlimited. 


Nonlinear  Inversion 

Robert  I.  Odom 
Applied  Physics  Laboratory 
University  of  Washington 
1013  NE  doth  Street 
Seattle,  WA  98105 

phone:  (206)  685-3788  fax:  (206)  543-6785  email:  odom@apl.washington.edu 

Award  Number:  N00014-08-1-0835 
http://staff.washington.edu/aganse/invresources 


LONG-TERM  GOALS 

The  long  term  goals  of  this  research  are  to  develop  practical  and  efficient  algorithms  for  application  to 
the  nonlinear  inversion  problems  encountered  in  ocean  acoustics.  Such  algorithms  would  be  used  for 
estimating  or  accounting  for  the  effects  of  the  environment  on  acoustic  propagation,  detection  and 
tracking  in  shallow  water. 

OBJECTIVES 

The  specific  objectives  of  this  research  are  to  adapt  a  specific  nonlinear  filter,  known  as  a  Daum  filter, 
for  acoustic  inversion  of  shallow  water  environmental  properties,  and  to  assess  the  performance  of  this 
nonlinear  filter  relative  to  local  linear  inversion  on  the  one  hand  and  global  methods,  e.g.  Monte  Carl 
methods  on  the  other  hand. 

APPROACH 

Many  inverse  problems  of  interest  in  ocean  acoustics  are  intrinsically  nonlinear,  e.g.  inverting 
measured  pressure  data  for  bottom  and  scattering  properties.  The  solution  to  the  nonlinear  inversion 
problem  is  usually  approached  in  one  of  two  ways.  The  first  way  is  to  assume  a  starting  model,  which 
one  hopes  is  near  to  the  true  model,  then  recursively  solve  a  linearized  version  of  the  inverse  problem 
for  corrections  to  the  starting  model  and  model  covariance.  The  advantage  of  this  approach  is  that  the 
numerical  implementation  of  the  solution  algorithm  is  relatively  straightforward  and  in  a  linear 
problem  the  statistical  properties  are  well  defined  and  will  remain  gaussian  if  they  start  out  gaussian. 
However  linearization  of  a  nonlinear  system  can  produce  biased  estimates  for  two  reasons:  1. 
Linearization  of  the  system  and/or  measurement  equations  may  not  be  a  good  approximation,  and  2. 
Nonlinear  systems  do  not  maintain  gaussian  statistics  as  they  evolve  even  if  they  are  initially  gaussian. 
Another  problem  with  linearizing  a  nonlinear  system  is  that  with  a  poor  starting  guess  the  solution 
algorithm  may  never  converge  to  the  true  answer.  If  the  starting  model  represents  a  point  near  a  local 
minimum  of  the  solution  space,  the  final  solution  will  be  trapped  in  that  local  minimum,  and  never 
converge  to  the  true  answer.  This  can  be  circumvented  by  using  Monte  Carlo  techniques  to  randomly 
sample  the  solution  space  for  starting  models. 
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The  other  class  of  solution  methods  attack  the  nonlinear  problem  directly  by  using  simulated 
annealing  or  genetic  algorithms.  The  disadvantage  of  these  directly  nonlinear  methods,  is  that  there  is 
no  way  to  conveniently  propagate  the  statistical  properties  of  the  solution  through  to  the  final  result. 
One  solution  to  this  problem  is  to  find  the  global  minimum  in  the  solution  space,  if  one  exists,  then 
linearize  about  the  solution  representing  the  global  minimum  and  do  a  statistical  analysis  about  that 
solution.  This  was  done  by  Potty  et  al.(2000),  who  employed  a  genetic  algorithm  followed  by  linear 
analysis  about  the  solution  determined  by  the  genetic  algorithm. 

The  recursive  algorithms  commonly  employed  for  the  estimation  of  the  model  and  covariance  relative 
to  some  initial  starting  values  bear  a  strong  resemblance  to  Kalman  filters,  which  are  commonly 
employed  in  target  tracking  algorithms.  The  original  Kalman  filter  was  derived  for  strictly  linear 
systems.  However,  the  Extended  Kalman  Filter  can  be  applied  to  systems  which  are  weakly  nonlinear. 
In  the  late  1980s  Frederick  Daum,  a  mathematician  working  at  Raytheon  Corporation,  developed  a 
fully  nonlinear  formulation  to  the  filtering  problem  for  target  tracking  (Daum,  1985,  1986,  1987).  His 
theory  is  elegant,  but  impractical  from  an  implementation  point  of  view.  Sometime  later  Schmidt 
(Schmidt,  1993)  succeeded  in  deriving  an  approximate  algorithm  based  on  Daum's  original  theory, 
and  developed  a  successful  numerical  implementation  of  a  nonlinear  filter  that  was  a  significant 
improvement  to  the  Kalman  and  Extended  Kalman  filters  for  the  type  of  tracking  problem  Schmidt 
was  interested  in. 

Filter  type  algorithms  are  ideally  suited  to  inverse  problems  with  time  dependent  oceanography  or 
range  dependence.  We  do  not  anticipate  attempting  to  include  time  dependent  oceanography  at  this 
time,  but  we  would  like  to  look  at  the  issue  of  range  dependent  inversion.  The  idea  would  be  to 
sequentially  update  parameter  estimates  as  a  function  of  range.  Also  note  that  any  inversion  algorithm 
can  be  cast  into  a  filter  like  algorithm  by  supplying  the  data  sequentially  and  updating  the  model 
parameter  estimates  sequentially  as  data  is  added  to  the  problem,  or  a  smoother  by  considering  the 
complete  data  set,  and  working  both  forwards  and  backwards  through  the  data  set.  In  the  end, 
probably  the  best  formulation  to  use  for  a  given  inverse  problem  depends  on  the  noise  statistics.  This 
is  also  something  we  propose  to  investigate. 

Finear  inverse  problems  admit  the  construction  of  both  data  and  model  resolution  matrices.  These 
resolution  matrices  can  be  used  as  metrics  with  which  to  estimate  model  uniqueness  and  data 
predictability.  We  will  be  able  to  construct  resolution  matrices  for  the  nonlinear  problem  and  compare 
them  with  their  fully  linear  equivalents. 

Quantification  of  the  resolution  of  an  inversion  can  be  used  for  experimental  design.  While  the 
resolution  of  a  linear  problem  is  well  defined,  and  described  in  basic  texts  such  as  Menke  (1983),  it  is 
less  so  for  a  nonlinear  problem.  One  of  the  objectives  of  this  research  is  the  quantification  of  the 
resolution  for  nonlinear  problems.  The  resolution  for  the  nonlinear  problem  can  be  defined  formally: 
Given  a  true  but  unknown  model  mtrue  such  that 

d  =  g(mtrue)  ->  mtrue  =  g‘^(d) 

where  d  is  the  data  vector,  and  g  is  the  operator  connecting  the  data  to  the  model,  e.g.  the  wave 
equation,  how  close  is  a  particular  estimate  mest  to  mtrue  ? 
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mtrue  =  gesf^(d)  =  gesf^[g(mtrue)] 
mest  =  r(mtrue) 

where  r(.)=  gest-'[g(.)]  is  the  model  resolution  operator.  The  nonlinear  model  resolution  r(.) 
operator  ean  be  eomputed  iteratively  from  the  Neumann  series  representation  for  gest  with  assumption 
that  both  the  data  functional  and  the  model  perturbation  functional  possess  regular  perturbation 
expansions.  (An  example  of  a  problem  which  possesses  a  model  perturbation  functional  with  a  regular 
perturbation  series  is  normal  mode  acoustic  propagation  with  “slow  enough”  perturbations  such  that 
the  modes  adjust  adiabatically  to  the  perturbations,  and  the  mode  eigenvalues  are  “far”  from  cut-off.) 

d  =  g(£m)  =  £G("')m  +  £2G(2)mm  +  £3G(^)mmm  +  ... 

£m  =  £l("')d  +  £2|(2)dd  +  £3|(3)ddd  +  ... 

where  the  l(i)  are  the  expansion  operators  of  the  g""'  Substitute  data  into  the  model  expansion  and 
order  by  £: 


£"':  m  =  l("')G("')m 
£2;  0  =  (lO)G(2)  +  l(2)G(1)G(1))mm 
£3;  0  =  ... 

The  remarkable  thing  to  note  about  these  expansions  is  that  the  nonlinear  components  of  the  model  in 

the  data  do  not  contribute  to  the  reconstruction.  (Snieder,  1990).  The  term  |(1)g(2)  is  a  linear 
inversion  of  the  component  of  the  data,  that  has  a  quadratic  dependence  on  the  data.  If  the  nonlinear 

inversion  is  to  reproduce  the  model  m  exactly,  the  )g(2)  term  must  be  canceled  by  the 

|(2)g("I  )g("'  )  term.  We  can  now  define  the  nonlinear  resolution.  For  the  estimated  model  mest  we 
have 

mest  =  l('')G('')m  +  (|('I)g(2)  +  |(2)F('')F(''))mm  +  ... 
mest  =  R("')mtrue  +  R(2)mtruemtrue  +  ... 


where  rC')  tells  us  how  much  smearing  there  is  in  the  map  between  m  and  mest,  and  R(2)  tells  us 
how  much  spurious  nonlinear  mapping  from  the  true  model  there  is  to  the  model  functional  {Snieder, 
1991). 
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WORK  COMPLETED 


This  year  we  focused  on  two  main  topic  areas  within  the  subject  of  continuous  geophysical  inversion 
of  the  ocean  bottom:  calculation  of  derivatives,  and  the  dependence  of  the  inverted  ocean  bottom’s 
resolution  upon  the  nonlinear  problem  formulation.  Both  are  briefly  summarized  below.  To 
adequately  model  the  geoacoustic  properties  of  the  ocean  bottom  we  use  the  elastic  wave  propagation 
paradigm.  While  some  investigators  solve  for  just  a  handful  of  parameters  describing  the  ocean 
bottom  properties,  we  instead  concentrate  on  the  continuous  inverse  problem,  solving  for  the  ocean 
bottom  properties  as  a  (possibly  piece-wise)  continuum.  Finite  information  in  the  measured  data 
certainly  limits  the  resolution  at  which  the  continuous  properties  can  be  solved  for,  but  one  can 
quantify  and  compare  this  resolution  as  well  as  the  associated  uncertainty  of  the  solution,  without 
inadvertently  imposing  false  structure  on  that  solution  (for  example  by  requiring  a  few  layers  when  the 
regional  structure  is  not  conclusively  known).  The  data  formulations  we  use  keep  the  problem  in  the 
weakly  nonlinear  domain  in  order  to  allow  a  linearization  approach  to  solving  the  inverse  problem. 
However,  derivatives  of  the  data  with  respect  to  the  bottom  properties  are  required  for  this  technique. 

Two  approaches  to  derivative  computation  were  explored.  The  first  was  analytical  formulation  based 
on  adjoint  methods.  The  starting  point  for  the  analytical  formulations  was  Tarantola’s  expressions  for 
these  derivatives  for  the  elastic  seismic  problem,  which  were  formulated  around  the  tensorial 
displacement  Green’s  function  of  the  elastic  problem.  However,  the  Green’s  function  calculated  by 
our  ocean  acoustic  codes  (still  with  an  elastic  ocean  bottom)  is  the  pressure  Green’s  function,  since  the 
source  and  receiver  are  in  the  water.  While  the  expressions  for  derivatives  of  just  the  P-wave  portion 
of  the  wave  propagation,  with  respect  to  the  P-wave  velocity  of  the  ocean  bottom,  could  be  shown  to 
match  previously  published  derivative  expressions  for  the  much  simpler  acoustic-only  case, 
reformulation  and  implementation  of  the  rest  of  the  derivative  expressions  for  the  full  elastic  problem 
is  work  in  progress. 

The  second  approach  to  derivative  computations,  which  we  finally  used  for  all  the  inversion  work,  is 
computing  first  finite  differences.  While  conceptually  the  simplest  approach,  it  is  very  computationally 
intensive  and  requires  analysis  in  the  choice  of  parameter  step  sizes.  To  make  the  method  feasible  in 
terms  of  computation  time,  distributed  computing  software  was  written  to  spread  the  computations  out 
over  many  processors  on  a  number  of  computers  around  our  lab  (varying  from  32-50  processors). 

RESULTS 

One  way  the  step  size  was  explored  by  plotting  (see  Figure  1)  the  difference  between  the  first  finite 
difference  values  computed  at  one  step  size  h  and  those  at  a  nearby  step  size  2h.  When  viewed  across 
orders  of  magnitude  of  step  sizes,  these  two  step  sizes  are  quite  close  and  ideally  should  yield  little 
difference  in  result.  However,  if  the  step  size  is  too  large,  one  increases  truncation  error  in  the  Taylor 
series  from  which  the  expression  for  first  finite  differences  is  derived.  If  the  step  size  is  too  small,  one 
increases  the  conditioning  error  caused  by  inherent  errors  in  the  G(m)  and  G(m-l-h)  terms  such  that  if 
their  values  are  very  close  then  their  difference  has  a  large  relative  error.  As  seen  in  both  plots  in 
Figure  1  the  optimal  step  size  lays  in  between.  An  example  of  the  form  that  these  derivatives  take  is 
seen  in  the  plot  in  Figure  2. 

Our  motivation  in  our  resolution  studies  is  to  explore  the  resolution  differences  for  different  data  and 
model  formulations  and  geometries  in  the  geoacoustic  inverse  problem.  Data  formulations  that  we  are 
comparing  include  travel-times,  tau-p  timeseries,  and  an  intensity  envelope  function;  model 
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formulations  (of  ocean  bottom  properties)  inelude  veloeities,  slownesses,  and  impedanees.  Geometries 
include  vertical  array  and  horizontal  array  configurations  at  different  source-receiver  ranges. 

Some  examples  of  these  quantifieations  are  shown  in  Figures  3  and  4,  taken  from  Ganse  and  Odom 
(2008).  Figure  3  shows  two  resolution  matriees  for  the  same  problem  exeept  for  differing  data 
formulations.  The  resolution  matrix  is  the  first-order  deseription  of  the  resolution  limits  of  the  inverse 
problem  solution  -  only  weighted  averages  of  profiles  ean  be  estimated  in  the  inversion,  and  these 
matriees  eontain  the  weights.  An  exaetly  diagonal  resolution  matrix  would  represent  perfeet 
resolution,  and  spread  from  the  diagonal  indieates  broadening  resolution,  or  smearing.  For  the  result 
on  the  left,  the  data  was  formulated  as  the  intensity  envelope  of  the  amplitude  time  series;  for  the  result 
on  the  right  it  is  formulated  as  a  tau-p  time  series  for  a  given  range  of  slownesses  (2.5e-4  to  8.0e-4 
s/m).  In  this  partieular  example,  the  tau-p  result  on  the  right  has  less  resolution  beeause  its  information 
was  windowed  in  slowness.  Figure  4  shows  geometrie  effeets  of  resolution  changes,  in  a  problem  with 
a  horizontal  array  near  the  oeean  surfaee  at  varying  ranges  from  the  near-surface  souree.  There  is  a  lot 
of  information  in  the  resolution  matrix;  using  the  traee  of  the  matrix  as  a  summary  value  we  eompute 
many  resolution  matriees,  one  at  eaeh  blaek  dot  on  the  plot.  One  ean  see  variation  and  maxima  in  the 
resolution  as  the  range  is  changed  in  the  problem. 


1  St  finite  dift  step  size  analysis  for  fwdprob-taup 
(parameters  are  d1  =1;  d2=40000;  m=1 0) 
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finite  diff  stepsize  h 
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finite  diff  stepsize  h 


1  St  finite  diff  step  size  analysis  for  f>vdprob~envelope 
(parameters  are  dl  =  10000;  d2=20000;  m=40) 


Figure  1:  Investigating  the  balance  between  truncation  error  and  conditioning  error  in  the 
computation  of first  forward  differences  of  the  ocean  bottom  geoacoustic  problem.  The  problem  has 
many  ocean  bottom  parameters  and  many  values  comprising  the  predicted  data,  so  any  given  plot 
can  only  show  a  snapshot  into  that  large  space  of  derivatives.  The  plot  on  the  left  here  analyzes  the 
derivatives  of  the  first  4e4  points  in  the  tau-p  time  series  with  respect  to  the  Itf'  ocean  bottom 
parameter  (p-wave  velocity  approx.  42m  below  ocean  bottom  interface).  The  plot  on  the  right  here 
analyzes  the  derivatives  of  le4 points  in  the  amplitude  time  series  envelope  with  respect  to  the  4(f'' 
ocean  bottom  parameter  (p-wave  velocity  approx.  166m  below  ocean  bottom  interface).  Note  the 
“happy  medium”  in  the  step  sizes  in  the  middle  of  each  plot  where  the  derivative  values  converge. 
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□  50  10D  15D 

depth  params  (2m  increments  to  300m  depth) 


Figure  2.  The  derivatives  as  approximated  by  first  finite  differences  of  the  amplitude  timeseries 
envelope  function  with  respect  to  p-wave  velocity  over  depth  in  the  ocean  bottom.  (Ocean  depth  = 
200m,  source  depth  =  10m,  receiver  depth  =  10m,  source-receiver  range  =  1km,  source  pulse  center 
frequency  =  lOOHz.)  Note  the  oscillating  nature  of  the  derivatives  caused  by  interference  of  up-  and 
down-going  acoustic  waves  in  the  vertically  stratified  medium;  the  amplitude  of  this  oscillation  dies 

out  exponentially  with  depth. 
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Figure  3.  Resolution  matrices  via  showing  differences  in  resolution  of  the  same  geoacoustic  inverse 
problem  (source  depth=10m,  HLA  depth=10m,  40  receivers  every  100m  from  100m  to  4.0km, 
lOOHz)  but  with  different  formulations  of  the  measured  acoustic  data.  For  the  left  result  the  data 
was  formulated  as  the  intensity  envelope  of  the  amplitude  time  series;  the  right  result  is  based  on  the 
data  formulated  as  a  tau-p  time  series  for  a  given  range  of  slownesses  (2.5e-4  to  8.0e-4  s/m).  Both 
the  horizontal  and  vertical  axes  on  these  plots  are  the  indices  of profile  depth.  The  resolution  matrix 
is  the  first-order  description  of  the  resolution  limits  of  the  inverse  solution  -  only  weighted  averages 
of  profiles  can  be  estimated  in  the  inversion,  and  these  matrices  contain  the  weights.  An  exactly 
diagonal  resolution  matrix  would  represent  perfect  resolution,  and  spread  from  the  diagonal 
indicates  broadening  resolution,  or  smearing,  at  a  given  depth.  This  particular  tau-p  result  on  the 
right  is  more  smeared  because  its  information  was  windowed  in  slowness. 
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Resolution  variations  with  geometry.  HLA  exampte,  5  revs  w/  100m  spacing 
Xdogspace(-6,-0.5,12),  smallest  X  has  largest  trace 


Resolution  variations  with  geometry.  HLA  example.  5  revs  w/  100m  spacing 
>.:=logspace(-7,-1.5,12).  smallest  X  has  largest  trace 


Figure  4.  Evolution  of  traces  of  the  resolution  matrices  (a  proxy  for  information  content  of  the 
inversion  solution)  as  a  function  of  HLA  receiver  array  range  and  regularization  parameter.  Each 
point  on  the  plots  corresponds  to  an  inverse  result  using  a  5-element  horizontal  array  of  receivers  at 
some  range  -  in  the  left  plot  a  density  profile  was  inverted  for,  in  the  right  plot  a  P-wave  velocity  plot 
was  inverted  for.  The  variations  seen  in  this  plot  are  not  strong,  but  one  can  see  trends  and  maxima 

in  the  resolution  as  the  experiment  geometry  is  varied. 

IMPACT/APPLICATIONS 

A  nonlinear,  well  eharacterized  inversion  method  and  algorithm  will  have  applieation  to  environmental 
estimation  and  target  tracking.  A  practical  method  way  to  compute  the  resolution  for  a  nonlinear 
inversion  will  have  an  impact  on  the  characterization  of  uncertainty  and  uniqueness  of  environmental 
estimates  required  for  acoustic  propagation. 

RELATED  PROJECTS 

Our  research  is  directly  related  to  other  programs  studying  effects  of  uncertainty  in  the  environment, 
measurements,  and  models  on  acoustic  propagation,  and  target  detection  and  characterization. 
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